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i-Q ■ Abstract. Free cooling of granular materials is analyzed on the basis of a pseudo- 

Liouville operator. Exchange of translational and rotational energy requires surface 
roughness for spherical grains, but occurs for non-spherical grains, like needles, even 
if they are perfectly smooth. Based on the assumption of a homogeneous cooling 
state, we derive an approximate analytical theory. It predicts that cooling of both 
rough spheres and smooth needles proceeds in two stages: An exponentially fast 
C/3 ' decay to a state with stationary ratio of translational and rotational energy and 

a subsequent algebraic decay of the total energy. These results are confirmed by 
simulations for large systems of moderate density. For higher densities, we observe 
deviations from the homogeneous state as well as large-scale structures in the veloc- 
ity field. We study non-Gaussian distributions of the momenta perturbatively and 
observe a breakdown of the expansion for particular values of surface roughness 
and normal restitution. 



1 Introduction 



■^- ' The hard-sphere model has been a very useful reference system for our un- 

derstanding of classical liquids Q. As far as static correlations are concerned, 
an analytical expression for the pair correlation is available 0] and provides a 
good first approximation for particles interacting via smooth-potential func- 
p^ ■ tions. The hard-sphere model is even more important for the dynamics, be- 

f^i cause it allows for approximate analytical solutions, based on the Boltzmann 

equation and its generalization by Enskog to account for a finite particle di- 
ameter and pair correlations at contact H. The model has the additional 
advantage that it is particularly well suited for numerical simulations pi and 
in fact many of the important phenomena of dense liquids have been observed 

^ ■ first in simulations of hard spheres. Examples are the discovery of long-time 

Q tails H and two-dimensional solids 0. 

O ■ Not surprisingly the model has become very popular also in the context 

of granular media, which are characterized by inelastic collisions of their con- 
stituents. Focusing on the rapid-flow regime, where kinetic theory should 
apply, generalized Boltzmann- and Enskog equations have been formulated 

C$ • and solved approximately ra,@,H,n0[O| . The success of the Boltzmann-Enskog 

equation in classical fluids is based on the linearisation of the collision oper- 
ator around local equilibrium. The resulting linear hermitean operator can 
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then be treated by standard methods of functional analysis |14^3||l^ |. For 
inelastic systems no analog of the local equilibrium distribution is known. In 
many studies, including the present one, a homogeneity assumption is made, 
which is known to be unstable for dense and large enough system and long 
times [15]. Hence the analysis is restricted to small and intermediate densi- 
ties. Alternatively, one may restrict oneself to almost elastic collisions and 
expand around the elastic case. 

Kinetic theory of rough, inelastic, circular disks was first discussed by 
Jenkins and Richman 0. These authors introduced two temperatures, one 
for the translational and one for the rotational degrees of freedom, and stud- 
ied deviations from a two-temperature Maxwellian distribution, using Grad's 
moment expansion. Subsequently Lun and Savage PjlOfl extended the ap- 
proach to rough, inelastic spheres. A set of conservation equations and con- 
stitutive relations was derived from the Boltzmann equation, assuming small 
inelasticity and surface roughness. Goldshtein and Shapiro |llf discuss in 
detail the homogeneous cooling state of rough spheres. They determine the 
asymptotic ratio of rotational to translational energy as a function of surface 
roughness and coefficient of normal restitution. Hydrodynamic equations and 
constitutive relations are derived with help of the Enskog expansion. More 
recently, event-driven simulations of rough spheres have been performed by 
McNamara and Luding [Q. They investigate free cooling as a function of 
arbitrary surface roughness and normal restitution and compare their results 
to an approximate kinetic theory |fl7||l8[ . 

Most analytical and numerical studies of kinetic phenomena have concen- 
trated on spherical objects so farg. The question then arises, which of the 
results are specific to spherical objects and which are generic for inelastically 
colliding particles. A single collision of two arbitrarily shaped, but convex 
objects is quite difficult to describe analytically pH , set aside the problem of 
an ensemble of colliding grains. In this paper we have chosen the simplest non 
spherical objects, needles, which allow for an analytical, albeit approximate 
solution and large scale simulations Q . 

The paper is organized as follows. In SecH we introduce the time evolu- 
tion operator. For pedagogical reasons we first discuss smooth potentials and 
recall the formalism of a pseudo-Liouville operator for elastic, hard-core colli- 
sions. Subsequently the formalism is extended to inelastic, rough spheres and 
needles. The homogeneous cooling state is introduced in Sec. 0. We present 
results for both spheres and needles, assuming a Maxwellian distribution for 
linear and angular momenta. We show with simulations that for dense sys- 
tems of needles the assumption of homogeneity breaks down. Corrections to 
a Gaussian approximation are discussed in Sec. ^. Finally in Sec. |5| we sum- 
marize results and present conclusions. Some details of the calculation are 
delegated to appendices. 



1 Exceptions are computer simulations of polygonal particles 119] and cellular au- 
tomata models [|2CJ] 
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2 The Liouville operator 

We are interested in macroscopic properties of systems of many particles 
which are themselves meso- or macroscopic, i.e. behave according to the laws 
of classical mechanics as opposed to quantum mechanics. In addition, our 
systems are granular so energy is not conserved. This means that they can 
not be treated with Hamiltonian mechanics. We will present here a formal- 
ism based on the Liouville operator that enables us nevertheless to derive 
properties of the system under consideration. 

We consider two different models: The first is a system of spheres of 
diameter d and the second is one of (infinitely) thin rods or needles of length 
L. In order to keep the discussion as transparent as possible, the formalism of 
the (pseudo) Liouville operator will be demonstrated for Hamiltonian systems 
with smooth potentials first, for hard core potentials next, and finally for 
granular spheres and needles. It is interesting to note that both cases, spheres 
and needles, arc analytically tractable so that comparisons between different 
geometrical particle shapes are possible. 

2.1 Smooth potentials 

We consider a system of TV classical particles of mass to in a volume V, 
interacting through pairwise potentials. The system is characterized by its 
total energy 

N 2 
i— 1 i<j 

in terms of particle momenta pi and coordinates r, . The time evolution of an 
observable /(-T), which is a function of phase space variables r := {ri,pi}, 
but does not depend on time explicitly, is given in terms of the Poisson 
bracket by 

f t ={H,f}=:i£f. (2) 

This defines the Liouville operator C The time evolution of / can then for- 
mally be written in terms of C: f{t) = e /(0). 

We decompose the Liouville operator C = £o + Antcr hito a free-streaming 
part Cq and an operator Anton which accounts for interactions. The definition 
of the Poisson bracket, 

{HJ} = ^{^-§-^)> (3) 
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with Tij = Ti — Tj . 



2.2 Elastic hard-core interactions 

A pseudo-Liouville operator for hard-core collisions has been formulated by 
Ernst et al. J23] and has been applied by many groups [g4j to study the 
dynamic evolution of a gas of hard spheres. Collisions are instantaneous and 
characterized by collision rules. In a collision of two particles, numbered 1 
and 2, their pre-collisional velocities v i = p\/m and v 2 = p-z/m are changed 
instantaneously to their post-collisional values v[ and v' 2 according to 

v[ =vi- {v 12 ri 2 )r 12 
v 2 = v 2 + (v 12 ri 2 )r 12 . 

We have denoted the relative velocity by V12 = v\ — t>2, and r = r/\r\. The 
free-streaming part of the Liouville operator remains unchanged, whereas the 
part which accounts for interactions has to be modified because the potential 
is no longer differentiable in the limit of hard-core interactions. As a conse- 
quence, £ is no longer self adjoint as it is for systems with smooth potentials. 
This is why it is called a pseudo-Liouville operator for hard-core systems. For 
the same reason we will need two Liouville operators below, one for forward 
and one for backward time evolution. 

In order to construct the pseudo Liouville operator, we consider the 
change of a dynamical variable due to a collision of just two particles. What 
we need is an operator 77 that 

• gives the change of an observable through a collision when integrated 
over a short time interval containing the collision time (since the hard 
core interaction is non-differentiable, we have to resort to integrating over 
the collision instead of looking at the derivatives directly) , 

• only acts at the time of contact, 

• only acts when the particles are approaching but not when they are re- 
ceding. 

The second requirement can be satisfied by T^ <x 5(\r i2\ — d), the third 
one demands T[ ' ex 0(— ^|ri 2 |), where ©(•) is the usual Heaviside step 

(12) 

function. In order to satisfy the first point, we use an operator b, which is 



defined by its action on an observable / according to 

b¥ ) f(vi,v 2 )=f(v' 1 ,v 2 ), (6) 
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(12) 



i.e. it simply replaces all velocities according to eqs. (|5|). The operator T| 

(12) (12) 

should give the change induced by a collision, so that T^ oc b\ — 1. 
We collect the three terms and make sure to include a prefactor which is 
chosen such that the integration of an observable over a short time interval 
around the collision time yields the change of the observable, as induced by 
the collision rules (^). The complete expression for T\ is thus 



iTl 12) = 



dt lri2 



6(\r 12 \-d)0(-±\r 12 \)(b^-l). (7) 



Since the probability that three or more particles touch at precisely the same 
instant is zero, we only need to consider two-particle collisions and find for the 
time-evolution operator for the system of elastically colliding hard spheres: 

f(t) = e l(Co+c±)t f(0) for i^ (8) 

with 



^ ± = E^ U) = E 



-i 

dr 



d , \ (y) 



1<J l<] 

da) 



S(\r Jt \ - d )9[ T f t \r j i\)(b^'-l). (9) 



The negative time evolution is given by £_, and b_ is the operator that 
replaces post-collisional velocities by pre-collisional ones. 



Extension to rough spheres Hard-core models of elastically colliding 
spheres have been extended to include rotational degrees of freedom and 
surface roughness flq,H. Rotational degrees of freedom offer the possibility 
to describe molecules with internal degrees of freedom and surface roughness 
is needed to transfer energy from the translational degrees of freedom to the 
rotational ones. 

We only discuss the simplest case of identical spheres of mass in, moment 
of inertia / and diameter d. Translational motion is characterized by the 
center-of-mass velocities «j and rotational motion by the angular velocities 
uJi. Let the surface normal fi2 at the point of contact point from sphere 2 
to sphere 1. The important quantity to model the collision is the relative 
velocity of the point of contact: 

V = (vi - -wi x f 12 ) - (v 2 + -u>2 x r 12 ). (10) 

There are two contributions, firstly the center of mass velocity of each sphere, 
and secondly the contributions from the rotations of each sphere. The minus 
sign in the first parenthesis stems from the fact that the surface normal, as 
it was defined, points outwards for sphere 2 and inwards for sphere 1. 
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Now we can specify the collision rules. Primed variables always denote 
quantities immediately after the collision; unprimed variables denote pre- 
collisional quantities: 

r 12 x V = -e t r 12 x V. 

As we are still dealing with elastic spheres, energy conservation requires 
e t — +1, corresponding to perfectly rough spheres, where the tangential 
velocity component is completely reversed. Perfectly smooth spheres et = — 1 
are also compatible with energy conservation, but reduce to the above simple 
case of spheres without rotational degrees of freedom, because during collision 
the angular velocities remain unchanged. Later, we will also admit other 
values for et- 

Eqs. (O) form three linearly independent equations. In addition, total 
momentum is conserved, 

v[ + v' 2 = vi + v 2 , (12) 

and forces during a collision can only act at the point of contact. Therefore 
there is no torque with respect to this point and consequently we have con- 
served angular momentum (also with respect to the point of contact) for both 
particles involved: 

vn d 

— f 12 x (v[ - Vl) + I(u}[ - U>l) = 
A ( 13 ) 

-^-r 12 x (v' 2 - v 2 ) - I(u>' 2 - u) 2 ) = 0. 

Altogether we have 12 independent equations for 12 unknowns, namely the 
four vectors v[ and u)\ with three components each. Solving for these, we 
obtain: 

v'x =vi- i] t v 12 - (ry„ - T)t)(ri2Vi 2 )ri 2 - i]t-r 12 x (u\ + w 2 ) 

v' 2 =v 2 + rjtw + {r] n - ■q t ){ri 2 v\ 2 )r 12 + rj t -r 12 x (wi + u> 2 ) 

Wi=Wi + -rVtri2 x t>i2 H r 12 x (r 12 x ui + w 2 )) 

dq q 

i 2 )i( A .„ , . . 

w 2 = o; 2 + —rjtru x t>i 2 H r i2 x (r i2 x wi + w 2 )). 

The dimensionless constant q = 41 /(md 2 ) abbreviates a frequently appearing 
combination of factors. We have also introduced two parameters r\ n and rjt, 
because we anticipate the more general collision rules for the inelastic case. 
For elastically colliding spheres, we simply have r\ n — 1 and r\t = g/(l + q) 
for perfectly rough and r/t — for perfectly smooth spheres. 
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The pseudo-Liouville operator for elastically colliding rough spheres is 
still given by eq. (M) but the operator b! now replaces linear and angular 
velocities according to eqs. dl4). 



Extension to rough needles Elastic collisions of hard needles have been 
discussed by Frenkel et al. Pq| . It is straightforward to rephrase their re- 
sults in terms of a pseudo-Liouville operator [Q . The free streaming part of 
the Liouville operator is derived from the kinetic energy of the Hamiltonian 
according to the general rules of classical dynamics. Note however, that for 
needles, one of the moments of inertia is zero; this implies that the angular- 
momentum component along the corresponding axis, which points along the 
orientation of the needle, is also always zero. Therefore, rotations about this 
axis can be ignored, and ui has only two components, both perpendicular to 
the orientation of the needle. The center of mass coordinate of needle i will 
be denoted by j", and its orientation by the unit vector Mj. The moments of 
inertia perpendicular to Ui are equal due to symmetry and will be denoted 
by /. 

The formulation of the collision rules proceeds in close analogy to rough 
spheres. First we determine the conditions of contact. The unit vectors U\ 
and «2 span a plane E\ 2 with normal 



u I = 



Ml X tt2 
\U\ X li 2 I 



(15) 



We decompose ri2 = r\— r 2 into a component perpendicular rj~ 2 — (ri 2 u±)u± 
and parallel r\ 2 = (s 12 iti — s 2 iu 2 ) to Eyi (see fig. 111). The rods are in contact 
if r 12 it ± = and simultaneously | S12 1 < L/2 and |s 2 i I < L/2. 



5 12 M 1 




Fig. 1. Configuration of two needles projected in the plane spanned by the unit 
vectors u\ and 112 
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The relative velocity of the point of contact is given by 

V = 1>12 + Sl2Ul - S 2 \U2- (16) 

It is useful to introduce a set of normalized basis vectors 



Mi, ttj 1 = (u 2 - (uiu 2 )u 1 )/y / l - («iM 2 ) 2 , and u± (17) 



with ttx defined in eq. (|15|). Total momentum conservation is given by ( |12[ ) 
and conservation of angular momentum with respect to the contact point 
reads 

u>i = u)\ H — Mi x (t^ — tq) and u; 2 = uj 2 H > — M2 x (u 2 — W2). 

(18) 

Three additional equations follow from the change in the relative velocity of 
the contact point, which is modeled in close analogy to the case of rough 
spheres: 

V'u± = -Vu±, V'ui = -e t Vui, and V'u 2 = -e t Vu 2 . (19) 

Again, energy conservation implies e t — ±1, corresponding to either perfectly 
rough or perfectly smooth needles (see also eq. (|32|)). Solving for v\ and u>' i7 
we obtain after a lengthy calculation: 

v[ = vi + Av and v' 2 = v 2 — Av (20) 

and u}[, u}' 2 given by eq. rtlq). The change in velocity Av can be decomposed 
with respect to the basis defined above, Av = 71M1 + 72M^ + au±. The 
coefficient a is given by 

2 ? \ ~ 1 

1 + ^f + ^fj Vu ± , (21) 

while 71 and 72 satisfy the set of linear equations 

AB\f ll \_ l + e t (V Ul \ (22) 



BC)\ l2 ) 2 \Vuf 

with 



■2 



A=i+^(i-(u lU2 n 



msn 



S = -^i(MiM 2 )v/l-(MiM2) 2 , (23) 



.2 



C = 1 + — - + — — (uim 2 ) 
2/ 2/ V ^ 
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The Liouville operator for two needles must obey the same basic require- 
ments as for spheres. The only changes are in the condition for a collision to 
take place, 

zT| 12) ex 9{L/2 - \s 12 \)0(L/2 - \s 21 \)6(\ri 2 \ - + ), (24) 

and in the condition that the two particles are approaching, 

l rr ) oce(-||^ 2 |). (25) 

Collecting terms and choosing the correct prefactor gives the result 



iT[ 12) = 



d 
dt lr 



"1-5"^ 



x 9(L/2 - \s l2 \)0(L/2 - \s 21 \)5(\r± 2 \ - 0+)(^ 12) - 1). (26) 
The operator b^ + ' replaces all velocities according to eqs. (^8|) and (|20|). 



2.3 Inelastic collision 

The collision rules for rough spheres and needles are easily generalized to 
inelastic collisions. This will allow us to set up a formulation of the dynamics 
of inelastically colliding grains in terms of a pseudo-Liouville operator. 

Rough spheres Energy dissipation is modeled by normal and tangential 
restitution. The collision rules imply for the change in the relative velocity 
of the points of contact: 

r 12 V> = -e n r 12 V 
r 12 xV = -e t r 12 x V. 

The first of these equations describes the reduction of the normal-velocity 
component by a non-negative factor e„. This is the well-known normal resti- 
tution. The second equation tries to describe surface roughness and friction 
in that it imposes a reduction or even a reversal of the tangential velocity 
component. It is motivated by the picture of small "bumps" on the surface 
which become hooked when the surfaces are very close. For all — 1 < e* < +1 
dissipation is present. 

The change in energy is given by 



AE= -? 



1-e 



2 



4 



-(r 12 v 12 ) 



2 



1 — e? q ( ,„ .„ eL / n\ 2 1 ,„„n 

— - — — — \vvi - {r 12 v 12 )r 12 - -ri2 x (u>i +u; 2 ) J . (28) 
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With the parameter range < e„ < 1 and — 1 < e t < 1, energy is only lost 
and never gained in a single collision. 

The conservation laws for linear and angular momenta are unchanged, so 
we obtain the same set of equations for the post-collisional velocities as eqs. 
14|), with however different parameter values 

1 + e„ q e t + 1 

^ = —2- and m= l + q 2 ' (29) 



Later we will need the inversion of eqs. (14), i.e. for given post-collisional 
velocities we want to determine the pre-collisional ones. This is simply done 
by replacing e t by l/e t and e„ by l/e„ in eqs. (114). The pre-collisional veloc- 
ities obtained from post-collisional ones will in the following be denoted by 
v",V2,uj'{ and w^ 

Rough needles For hard needles we introduce normal and tangential resti- 
tution according to: 

V'u± = —e n Vu±, V'ui = — etVui, and V'u 2 = —e t Vu 2 . (30) 

The conservation laws for linear and angular momenta are the same as for 
the elastic case, so that one arrives at the same set of eqs. (p0[), the only 
change affecting the parameter 



a =-i±-( 1+ ^ + ^yV^ pi, 



The energy loss for needles is given by 

1 - ef fC(V Ul ) 2 - 2B(Vu 1 )(Vui) + A(Vui) 2 



AE= -? 



AC-B 2 



„2 - '' ' • ' 



m l_£»(, + ^i + ^i] {Vu _r. ,-;,, 



It can be checked with eqs. (E3J) that the first term is less than if and only 
if — 1 < et < 1. Obviously, the second term is also less than if < e„ < 1. 
Our method of modeling granular collisions of needles is therefore consistent 
with the constraint that energy may not be gained in a single collision. 

2.4 Time evolution of the distribution function 

We will be interested in ensemble averages of observables /(/") at a time t 
defined by: 



(/}(*) = J dr p(r-o)f(r-t) = J dr P (r ; t)f(r). 



(33) 
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Here p(-T; t) is the TV-particle distribution function at time t. The average can 
either be taken over the initial distribution p(r; 0) at time 0, the observable 
being propagated to time t, or equivalently over the distribution p(-T; t) at 
time t with the unchanged observable f(r). We write eq. (p3|) as 



(34) 



(/>(*) = dr P (r;0)e' Lt f(r) =■. dr U Lt P {r-Q) /(r), 



to define the time-evolution operator C which describes the time evolution of 
p |2q] . To determine C explicitly, we take the derivative of eq. (pj) at time 
t = for simplicity, 



9t(f)(t)\ t=0 = dr P (r;0)i£f(r) 



dr (d tP (r ; t)\ t=Q ) f(r) = / dr (*£p(r ; o)) f(r) 



(35) 



The time-evolution operator of the density due to free streaming, Lq, is 
easily calculated by partial integration and we get Cq = —£q. To find an 
expression for the time-evolution operator of the density due to collisions 

T + for spheres, we use eq. (p5[). Phase-space coordinates before collision 
are denoted by r, after collision by r' 



6i 12) r so that 



f dr P {r- o) ? r{ 12) 



an 



I 



dr P (r;0)S(\r 12 \ - d)0(--\r 12 



dt 



ri2 



(f(r')-f(r)). (36) 



In the first term on the right hand side we make a coordinate transformation 
to the variables after collision with Jacobian J :— |^p-|- We use the inverse 

operator of b. 



namely b_ F = T . Here the coordinates before collision 



in terms of the coordinates after collision are denoted by T" — r(r'). We 
^12^12 an d rewrite the first term 



note that 2§l r i2 



dr P {r- o)<5(|n 



d)0(--\r 12 \ 



dt 



ri2 



f(r') 



J dr'Jp(r";t)6(\r 12 \ -a)0(-< 2 f 12 )K 2 f 12 |/(r') (37) 

Next we rename J" by r and make use of v'^ m r nm = — —(v nm r nm ). This 
allows us to identify the time-evolution operator of the distribution function, 

T+ , by: 



iTf ) = <5(|r 12 | - d) 



dt^ 



9Ar 12 \)fbW-9(-±\r 12 \)). (8s) 
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It is common to rewrite eq. ( pq ) by multiplying it with J dcr8(cr — Vyi) 
so that we can replace r\i by <x in eq. (pq). In the second term the integral 
transformation <j — > — er is done and we integrate over \u\. We obtain in D 
dimensions 

iTf ] = d°- x f da (v 12 a) (—S(r 12 - da)b^ 2) - 6(r 12 + da) 

Jv 12 a->0 \ e " 

(39) 

Finally, we note that t = is not special since we have only chosen it for the 
sake of simplicity. Hence we have derived the time-evolution operator for the 
iV-particle distribution function p{T; t) which is given by the pseudo-Liouville 
equation 



d tP (r,t) = i -coin +E T + J) p^ 1 )- ( 4 °) 



Mil) 
i<j 



A similar procedure yields the time evolution operator for the distribution of 
needles. 



3 Homogeneous cooling state 

We are interested in the time evolution of a gas of freely cooling rough spheres 
or needles which is dominated by two-particle collisions, as discussed in the 
previous section. We aim at a description in terms of macroscopic quantities 
and focus on the decay in time of the average kinetic energy of translation 
and rotation, defined as 

(Etr) (t) = ^Y,[ dr P( r ^ V * =■ if TtrW ' 
i 

&<*) (*) = j^Y,J dr p( r ^ =■ ^f T ^(t) ■ 



(41) 



Here D tT and D mt denote the total number of translational and rotational 
degrees of freedom respectively. It is impossible to compute the above expec- 
tation values exactly and we have to resort to approximations. We assume 
that the A^-particle probability distribution p(r, t) is homogeneous in space 
and depends on time only via the average kinetic energy of translation and 
rotation: 

pacs(r;t)iic8~W(r 1 ,...,r N )p({v i ,u i }iT tt (t),T tot (t)) . (42) 

The function W(ri, . . . , Tat) gives zero weight to overlapping configurations 
and 1 otherwise. Needles have vanishing volume in configuration space, so 



Free cooling of particles 13 

that W = 1. We shall furthermore assume that p factors neglecting correla- 
tions of the velocities of different particles. In the simplest approximation we 
take p to be Gaussian in all its momentum variables 



p({Vi,u}i};T tI (t),T Iot (t)) ex exp 



N i E tr , E IC 



Tt r (t) Trot(i) 



(43) 



In the next section, we shall discuss non Gaussian distributions and shall 
compute corrections perturbatively. 

To determine the time dependence of Tt r (t) and T ro t(i) we take time 
derivatives of eqs. §H\) and use the identity f t (f) (t) = J dr{f t p(r, t))f(r) = 
J dr{iCp{r,t))f{r) = f dr P (r,t)iCf(r). Then p(r,t) is replaced by phcs(IV), 
resulting in 



d 2 f 2 

— T tI (t) = -—- / dr p HC s{r;t)iCE tl . = —— (i£E tI ) HCS and 

at JJtr J Mr 

d 2 f 2 

-r;T mt (t) = — — / dr p HCS (r;t)i£E Iot = — — (iCE r 
dt D mt J D mt 



(44) 



ot/HCS 



All that remains to be done are high dimensional phase-space integrals, the 
details of which are delegated to appendices |A| and [B|, for spheres and needles. 

3.1 Results for spheres 

After integration over phase space has been performed (see Appendix A for 
details), we find 

%^T tr (t) = (i£E tI ) ucs = -GAT*' 2 + GBT?/ 2 T rot , 



2 dt 

-Prot d 

2 dt 



— -T mt (t) = (i£E mt ) ws = GBT^J 2 - GCT t \ /2 T mt , (45) 



with the always positive constants A, B, C ', and G depending on space di- 
mensionality D. In two dimensions the constants in eqs. (Eq) are given by 

G = 4d^[^g(d), A = i-^ + f (1 - m ), 

V V m 4 2 

c -tH)- '«' 



A= 1 -^+ Vt {l-r lt ), 
C=3*(l-5*V (47) 
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2q 
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The pair correlation function at contact, g(d), is denned in the usual way G|. 
A detailed discussion of these results, and in particular the dependence of 
free cooling on e„ and et, can be found in fll8| . 

The Enskog value |14|,|28| of the collision frequency uje, i-e. the average 
number of collisions which a particle suffers per unit time in D dimensions is 
given by 



u, E :=S D ^ 9 (d)d D -'J^. (48) 

V V irm 

Sd is the surface of a unit sphere in D dimensions. Note that always uje oc 
GT 1 / 2 . We define dimensionless time r by dr = u E dt so that r counts the 
collisions that on average each particle has suffered until time t. In a simu- 
lation this would simply be done by counting the number of collisions. The 
functional dependence of the two temperatures on r is determined by 

-^-T tr = -aT tl . + 6T rot , (49) 

dr 

-—Trot = bT tr - cT rot (50) 

dr 



with properly defined a,b,c. Eq. (49) has a simple interpretation: In a given 
short interval At a number of At collisions occur. Due to these collisions 
translational energy decreases by an amount given by the first term, but 
there is also a gain term, reflecting that rotational energy is transfered to 
translational energy. The solution of eqs. ( |49| , [50|) can be written as 

T tr = Cl K+ exp(-A+t) + c 2 K_ exp(-A_£) , (51) 

T Iot = c\ exp(-A+t) + c 2 cxp(-A_i) , (52) 

1 



K ± = -{c~a±^{c-ay+AV) , (53) 

A± = i(c + aT V(c-a) 2 +4& 2 ) . (54) 

The constants c\ and c-i are determined by the initial conditions and A_ > 0, 
A + > and A_ > A + holds for all e*, e„. Hence for long times the ratio of 
7tr /Tot is determined by K + . 

We now assume that the ratio T tI /T TOt has reached its asymptotic value 
K + for some r > t or equivalently t > t and substitute T rot = T tI /K + into 
eq. p3) we obtain 

±T tI = -FT 3 / 2 . (55) 

The resulting equation is of the same functional form as for homogeneous 
cooling of smooth spheres, except for the coefficient F, which contains all the 
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dependence on system parameters. Its solution is given by 

T tr (io) 1 



Ttr 



[l + T tI (toy/ 2 (F/2)(t-t )} 2 {Ft/2f 



(56) 



Haff 's |27| law of homogeneous cooling. We have determined two time scales, 
first an exponentially fast decay (measuring time in collisions) towards a 
state where we find a constant ratio of translational and rotational energy. As 
long as dissipation is small, we approximate the Enskog-collision frequency 
for sufficiently short times by its initial value u)(t) ~ w(0) so that we find 
exponential behavior also in real time. The second stage of relaxation is 
characterized by a slow, algebraic decay of both energies, such that their ratio 
remains constant. These two time regimes are clearly seen in the numerical 
solution of eq. ( Eg ) for initial conditions T tr (0) = 1 and T ro t(0) = 0, i.e. 
a system prepared in an equilibrium state of perfectly smooth spheres. We 
show in fig. a) in a double logarithmic plot the time dependence of the total 
c 1 -J(T tr (t)+T rot (i)) and the ratio T tr (i)/T rot (t). Time is plotted in 

0.9 and e t = -0.8. 



energy E 



units of lGT t / (0). We have chosen e 




10 10 10 10 10 10 10 10 10' 



b) 




rt 



Fig. 2. a) Theoretical prediction (lines) for spheres for the total energy E — 
§(r tr (t) +T ro t(i)) and the ratio Tt T (t)/T TO t(t) versus time. Time is plotted in units 
of |GT tr (0). We have chosen e n = 0.9 and et = —0.8. The symbols represent data 
of a simulation of 1000 particles in a box of length 16 d. 

b) The same as a) for needles: The total kinetic energy E — 3/2T tr + T rot and 
Tt r /T r ot are plotted vs. time in units of 7n\/Ttr(0). The simulation data are from 
a system of 10000 needles in a box of length 24 L with e„ = 0.8. 



3.2 Results for needles 

In the case of needles we restrict ourselves to the case of perfectly smooth 
needles, i.e. et = — 1. After some lengthy algebra, presented in appendix pi 
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eq. (H) can be cast in the following form 

2T tr _ f 2 (1 + ^fcr 2 ) 1 / 2 



= - d 



7 „T t 3 / 2 (l + e n ) •/□ ' 1 + kr2 



11 r C\ -L T Iot l, r 2\3/2 

1 + e "'«^T^' <-) 



□ 



4T rot /■ . ^/cr^l + ^fcr 2 ) 1 ^ 

- ' a r- 



3 7 „T t 3 / 2 (l + en ) 7d 1 + fcr 2 

i±£!i /■ fcr 2 (l + ^fcr 2 )3/ 2 
+ 2 ] n dr (1 + kr 2 ) 2 ' ( j 

with 7„ = (2JVX 2 A /7r)/(3V v / m) and fc = (mL 2 )/(2I). The two dimensional 
integration extends over a square of unit length, centered at the origin. 

In fig. g b) we plot the numerical solution of eqs. ( p7J , p8|) for e n — 0.8 and 
k — 6 (k — 6 corresponds to a homogeneous mass distribution along the rod) 
as a function of time in units of 7„ y/TtJfi). In addition we have performed 
simulations of a system of 10000 needles, confined to a box of length 24 L. 
We show the total kinetic energy E = |T tr + T rot (in units of T tr (T — 0)) and 
the ratio T tr /T rot . Analytical theory and simulation are found to agree within 
a few percent over eight orders of magnitude in time. (T ro t(0) = has been 
chosen as initial condition) . For needles we observe an even clearer separation 
of time scales. The decay of T tT /T m t to a constant value K + happens on a 
time scale of order one. In this range of times the total kinetic energy E 
remains approximately constant (on a logarithmic scale) and decays like t~ 2 
only after translational and rotational energy have reached a constant ratio. 
We plug the ansatz K + T rot = T tT into eqs. (pTlpq) and recover Haff's law 
also for needles. To determine the constant K + we use K + T lot — T tI = 0, 
which yields an implicit equation for c. Equipartition holds for all values of 
e„ if k = (mL 2 )/(2I) is set to particular value (k* — 4.3607), given as the 
solution of 

(1 - e 2 ) / d 2 r l ~ 2k * r2 = . 
V nJ Ju Vl + k*r 2 

For k < k* we find T tI < T rot and for k > k* , T tI > T rot . Hence the distri- 
bution of mass along the rods determines the asymptotic ratio of rotational 
and translational energy, including equipartition as a special case. 



3.3 Breakdown of homogeneity in dense systems of needles 

It is well known that dense and large systems of inelastically colliding spheres 
exhibit clustering so that the assumption of homogeneity breaks down and 



Free cooling of particles 



17 



deviations from Haff's law of homogeneous cooling are observed |l3,pG|. To 
investigate inhomogeneities for dense systems of needles we measure hydro- 
dynamic quantities, i.e we define local variables as the density field, the trans- 
lational and rotational flow field and the local rotational and translational 
kinetic energy. In order to take local averages over small volumes, we divide 
the simulation box into cells whose sizes are small compared to the box size 
but large enough to give a reasonable statistics. We choose cell size, such 
that on average about 25 needles are in each cell. For each cell indexed by 
a we compute the number density p a := pJ— X^eceli 1 = (1)"' the trans- 
lational energy per particle p a E^ — (J^vf) a and the hydrodynamic tem- 
perature T^ = E^ — mU 2 a /2 defined by fluctuations around the flow field 
PcJJa = (v) a - The corresponding observables of the rotational degrees of 
freedom are the rotational energy per particle ££ ot the hydrodynamic rota- 
tional temperature T£ ot and the rotational flow field fi a . 

To check for spatial clustering, we compare the statistics of fluctuations 
of the local density, velocity and translational energy for elastic and inelastic 
systems. As an example we show in fig. a the histogram of the deviation of 
the local density Sp a — p a jn — 1. We performed simulations of a dense and 
large system of 20000 needles confined to a volume with linear dimension 
12L and e n = 0.9. The initial distribution is uniform, corresponding to the 
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Fig. 3. Histogram of density fluctuations in the initial state and after 560 Collisions 
per particle. It is obvious that regions with high density have developed. 



equilibrium state of an elastic system. As the system develops in time with 
particles colliding inelastically, we observe that the distribution broadens, a 
clear indication that regions of large density have developed. Histograms of 
the local translational and rotational energies look very similar. 

Inelastic hard spheres without surface roughness tend to move more and 
more parallel so that large scale structures in the velocity field develop. 
In such a state most of the kinetic energy is to be found in the energy 
of the flow field, whereas the energy of the fluctuations around the flow 
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field is small. A quantitative measure for this effect p9| is the ratio of 
the total energy of the flow to the total internal energy of fluctuations: 
Str := (Y,aTPa U a)/(I2 a Pa T a) and tne analogous quantity S mt for the 
rotational degrees of freedom. In fig. we show St r and S TO t as a function of 
time, measured in collisions per particle. We observe an increase of St r by a 
factor of 50, whereas S rot increases only by about 50 %. Hence the large scale 
structures in the flow field are much more pronounced for the translational 
velocity. In the fig. || we show the flow field after 600 collisions per particle. 



0.08 
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Fig. 4. Ratio Str (Smt) of the local macroscopic energy to the local temperature 
for the translational (rotational) degrees of freedom as a function of the number of 
collisions per particle. 



We observe two shear bands (note the periodic boundary conditions) in which 
the flow field is to a large degree aligned. In periodic boundary conditions 
stable shear bands have to be aligned with the walls of the the box. 




Fig. 5. Flow field after 600 collisions per particle 
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How does the organization of the flow field influence the decay of the av- 
erage energy in the system? Brito and Ernst [B(j have suggested a generalized 
Haff 's law to describe the time dependence of the kinetic energy of smooth 
inelastically colliding spheres even in the non-homogeneous state. They found 
that in the late state where one finds a well developed flow field the energy 



decays like t~ d ' 2 in D dimensions. As in section 3.1 r is the average number 
of collisions suffered by a particle within time t. In fig. M we compare the data 
of the simulation with the solution of eqs. (pT^pq) and in the inset we plot T tr 
as a function of r and compare it to r -3 ' 2 . We can not confirm a r -3 ' 2 law, 
but by inspection of fig. |fj we see that the range of correlations are already 
of the order of the system size, so that finite size effects - not taken into 
account in the theory of Brito and Ernst - may be dominating. To simulate 
larger systems and longer runs has not been possible because simulations of 
dense systems are rather time consuming G2] . 
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Fig. 6. Data of simulations for translational and rotational temperature as a func- 
tion of 

(PIP!) 



time (units of 7nv/Tt r (0)) are compared to the numerical solution of eq. 
and to t -3 ' 2 . The inset shows T tr as a function of r and r -3 ' 2 . 



4 Non-Gaussian Distribution 



In this section we keep the assumption of homogeneity and factorization of 
the ./V-particle distribution function, but go beyond the approximation of a 
purely Gaussian state. Initially the system is prepared in a Gaussian state, so 
that deviations from the Gaussian should be small for short times and per- 
turbation theory can be used to check the range of validity of the Gaussian 
approximation. We expand the one particle distribution function in general- 
ized Laguerre polynomials (for a definition see pW) around the Gaussian with 
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time dependent variances. We define an average velocity vo(t) = \j2T tr {t)/m 
and u>o(t) — y // 2T I0 t(t)/I and scale linear velocities by Vo(t) and angular ve- 
locities by u>o(t). The general ansatz for the N-particle distribution function 
of the homogeneous cooling state then reads 

N 

p(r,t) ~ W(r 1 ,...r N )Y\_Pi(vi,Ui,t) , and 

1 
Pi(Vi,(ji}i,t) = — -— cxp 




(59) 



n,m=0 



We have introduced the abbreviations a — Z?tr/2 — 1 and (3 — D ro t/2 — 1. 
The average linear and angular velocities, Vo(t) and uio(t), are time dependent 
and so are the coefheients a n ,m(t) of the double expansion. At time t = 
the system is equilibrated with temperature T so that ^Vq = %T and 



jw 2 = ^f-T and hence a n , m (t) = 



20 2 

The factor i?(£) follows from the proper normalization, J dviduJiPi = 1, 

Z(t) = t^'u^V^V^Vo , (60) 

and we require that t>o(£) and u)o(t) be determined by the conditions 

J drv\p{r,t) = ^ 2 (t) and J drujp(r 7 t) = ^u; 2 (t) . (61) 

The orthogonality relations of the Laguerre polynomials imply ai,o(i) = 
«o,i(i) = f or a ll times t and 

<wW = {n l a) ^ J dr P (r,t)K ((^) 2 ) l£ ((^) 2 ) . (62) 

The binomial coefficients are denoted by (?) and we choose ao.o = 1- 

Taking the time derivative of eqs. ( |6l| , |62| ) one gets the full time dependence 
of the homogeneous cooling state given by the time dependence of all its 
momenta. Taking time derivatives of the right hand side of eq. (p2) one has 
to take into account the time dependence of p(r,t), which is determined by 

C as well as the time dependence of L" ((f 1 ) 2 ) Lf^ n ((^ i ) 2 ) via Vo(t) and 

uJo{t), which follows from eq. (|6l|) . 

Assuming that all a n , m are stationary in time and that vq/luq = p is 
constant we get an infinitely large, nonlinear system of equations. To make 
further progress we truncate the expansion in eq. ( J59J ) and take into account 
only a n ,m for n + m < 2. We also neglect in the system of equations products 
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of different a n _ m , which we assume to be of higher order. We show results 
for ao,2 in fig. ffl for fixed e n = 0.9 as a function of et- Deviations from 
the Gaussian vanish for perfectly smooth spheres and are found to increase 
dramatically for et — ► —0.9. Deviations from the Gaussian distribution are 
also small for perfectly rough spheres which is unexpected, because rotational 
degrees of freedom are coupled to translational ones and e n — 0.9. In fact 
deviations stay small for a broad range of values of e t >^ —0.75. We don't 
consider it meaningful to plot the theoretical result, once a divergence of ao,2 
has occurred. We measured ao,2 in simulations of small systems. Thereby we 
avoid clustering but have to bear with poor statistics. The simulations confirm 
the increase of ao,2 around et = —0.7 in agreement with the perturbation 
expansion. 

Goldshtcin and Shapiro flTl| ] propose a similar set of momentum equa- 
tions but they solve it only to lowest order, resulting therefore in the same 
asymptotic ratio /i as given in eq. (p3|). 




Fig. 7. Coefficient ao,2 for e n = 0.9 as a function of et- Theory (straight line) and 
simulations (circles) are compared. 



5 Conclusion 



Two simple models of granular particles with rotational degrees of freedom 
are discussed: Rough spheres or discs and, as an example for non-spherical 
particles, needles. We focus on the simplest collision rules, which allow for a 
transfer of translational energy to rotational degrees of freedom. For spheres 
this is achieved by tangential restitution (in addition to normal restitution), 
for needles normal restitution is sufficient. We show that the time evolution 
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can be formulated in terms of a pseudo Liouville operator, thereby gener- 
alizing previous work on elastic collisions to inelastic ones. The presented 
formalism is general enough to include more realistic collision rules, for ex- 
ample Coulomb friction for small angles of impact and tangential restitution 
for large angles. Work along those lines can be found in [B2| . 

The computation of non-equilibrium expectation values, like e.g. the re- 
laxation of kinetic energy, require approximations. These are formulated for 
the TV-particle distribution function, which we assume to be homogeneous 
in space and depend on time only via the average translational and rota- 
tional energy, T tr and T rot . The distribution of linear and angular momenta 
is expanded in Laguerre polynomials around a Gaussian state with time de- 
pendent width, Ttr and T rot . The zeroth order approximation, i.e. a pure 
Gaussian, leads to two coupled differential equations for the two tempera- 
tures. In both systems, spheres and needles, the relaxation of translational 
and rotational kinetic energy is characterized by two time scales: (1) An ex- 
ponentially fast decay towards a state with constant ratio of translational to 
rotational energy and (2) an algebraically slow decay of the whole energy, 
such that the above ratio keeps constant in time. The theoretically predicted 
cooling dynamics is supported by computer simulations of systems of small 
or moderate density, where no shearing or cluster instability is observed and 
the system remains homogeneous ]18| , [22| 

To study deviations from the Gaussian state, we restrict ourselves to 
rough spheres and truncate the expansion in Laguerre polynomials, keeping 
the first three terms (the first order term for smooth spheres has already been 
computed in |33|). This perturbative approach is shown to break down for 
certain values of e* and e ra , where deviations from the Gaussian are shown 
to diverge. These results are confirmed by simulations. We indicate, how a 
more general expansion with time dependent coefficients can be achieved. 
For totally smooth spheres this expansion has been performed up to fifth 
order |j4j . It predicts an exponentially fast decay of the coefficients to their 
stationary values in agreement with simulations and direct solutions of the 
Boltzmann equation |35| . 

For needles we observe and investigate the breakdown of homogeneity in 
simulations of dense systems, where the inter-particle spacing is smaller than 
the length of the needles. Large-scale structures in the translational velocity 
field are seen to develop. Furthermore the density does not remain homoge- 
neous but clusters form and dissolve again. These effects lead to deviations 
from the solution of the homogeneous cooling state on the longest times scales 
and a third stage of cooling is found. It is characterized by an even slower de- 
cay of the kinetic energy, most of the energy being stored in the macroscopic 
velocity field. 

We plan to derive generalized hydrodynamic equations for grains with 
rotational degrees of freedom and in particular hard rods. Such a set of hy- 
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drodynamic equations could serve as a starting point for a stability analysis, 
similar to the work of Brito and Ernst |30| for smooth spheres. 

Acknowledgement— This work has been supported by the DFG through 
SFB 345 and Grants Zi209/5-l and Zi209/6-l. 



A Calculations for spheres 



In this appendix we explain, as an example, the main steps to calculate 
(i£2£tr)HCS °f ec l- (0) m ^D. We define the configuration integral 

N 

Q N := j[[driW(ri,...,r N ) . (63) 

The proper normalized JV-particle distribution function for the HCS-state 
reads 

1 / \ N / T \ N / 2 



(64) 



pucs(r;t) = - w W(r u ...,r N ) [ „_ m ^ ) ( _^ M ) 



exp 



27rT tT (t)J \2TrT mt (t) 

N / T 

El m 2 1 2 



The angular velocity is a scalar in two dimensions, but a vector in more than 
two dimensions. Free streaming does not change the energy, so we have to 
take into account only the collision operator £ + and compute 



(iC + E tI ) ncs = W [dr Pws (r ; t^±Y,j v t = 
ijtj J fc=i 

The binary collision operator 7^ gives a contribution only, if either k = i 
or if k = j. Next, we introduce two ^-functions, 

{iC+E tI ) ws = ^J2J dr J dR 1 dR 2 6(R 1 - r l )5{R 2 - r 3 ) 

PHCsir-^iT^ivl + v*) , (66) 
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which allows us to replace r^ by i?i and 7*j by R 2 in Tp . We define the 
pair correlation function g(]R\ — -R2I) by 



^(l-Ri--R 2 | 



1 if 

nT,qn I[dr k W(r 1 ,...,r N )5(R 1 -r i )5(R 2 -r j ). (67) 

Eq. ( p7| ) is used to rewrite eq. (|66j) in terms of the pair correlation function. 
Integration over all velocities and angular velocities with index k and i ^ k ^= 
j gives 1 due to normalization. We get 



N , .2 



flf(r) |t>ia • r\ {-v 12 ■ r) 5 (|r| - d) 4£ tr . (68) 

The loss of translational energy of two colliding particles is denoted by AE tI . 
We use the abbreviation R\ — R2 = r = rr and neglect non contributing 
terms linear in fl so that AE tr is given by 

AE tI = ^[2 Vt (vt ~ 1)(«? 2 - («ia • rf) - 

(1/2)(1 - e»)(t;i2 • r) 2 + (1/2)^^^ + uo 2 f] . (69) 
To perform the remaining integrations we substitute 

f2 = — =(wi +ui 2 ), u> = — =(wi-w 2 ), (70) 

r = i2i-i? 2 , R = Ri. (72) 

The Jacobian determinant for the above transformation is 1. Integration over 
w, V and R all give the value 1 due to normalization. We are left with 

N m ( 21 \ 1/2 f 

exp {~^m - 2^iy ) ff(r) '" ■ f ' e ( -" • f } * (|p| - d) 



- [2 Vt ( Vt - 1)(« 2 - (v ■ f) 2 ) - (1/2)(1 - e 2 n )(v ■ rf + (l/2) V ?d 2 tf 
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The integration over \r\ yields dg(d). Choosing e.g. r to point along the x- 
axis, the integrals over linear and angular velocities can easily be done as 
moments of a Gaussian distribution. The result is independent of f , so that 
the integration over r gives 2tt. Finally we obtain the result of eq. (ffq). 

B Calculations for needles 

In this appendix we present some of the detailed calculations for needles. 
As a first step we express the orientation of the rods in spherical coordi- 
nates Ui = (sin(^) cos(0i),sin(6'i) sin(0i), cos(6*i)). The canonical momenta 
(translational and rotational) are then given by 

p t = mvi , pg z = I9i , p^ = I<j>i sin 2 9 . (73) 

In the following calculation it will be necessary to express iii in terms of 
canonical momenta 

eg i and e ( f >i are orthogonal unit vectors in 9i and 4>i direction. The kinetic 
energies per particle are then given by 

1 N 1 1^1 1 

e « = nY,^pI e «* = nY,vpI + ^-^pI- ( 75 ) 

i=l i=l ' 

We want to calculate non-equilibrium expectation values with the normalized 
probability distribution given in eq. (^3J). We consider again as an example 
the translation energy per particle Et r . 

IT P > l l l l 



V N (4tt)^ (27rMT trans )3JV/2 (27r/T rot )^ 
1 f N 

2 x! / n dr i d ^ d0 i dp i dpe 3 dp ^ 



exp{-NE tT /T tI (t) - NE Iot /T Iot (t)]iTl nm) E tI . (76) 



operator T+ gives a contribution only, if either i = n or if i = m. We can 



Similar to the calculation for the spheres we see that the binary collision 
operator f^ nm ' gives a contribution only, if eit 
sum over N(N — 1) identical integrals and get 

2 

N 2v' (If (2 Jr te )3 (2nITrotr /jg^ ^ ^ ** *»> ** 

exp[-^ a /rtr(t) - ^rot/r rot (t)] 



d 

dt 

0(X/2 - | Sl2 |)0(i/2 - |a a i|W|r£,| - 0+)Z\< . (77) 



d | x, 
dt I 12 ' 



e(~K 2 ) 
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E\^ (£rot ) i s the sum of the translational (rotational) kinetic energy of parti- 
cle 1 and 2 and with AE\^ we denote the change of the translational kinetic 
energy of particle 1 and 2 in a collision: 

AE% = (Pl ' P2) ^ P + *£- , (78) 

"mm 

Ap = _I+£z_ ^^(V ■ u ± )u ± . (79) 

m ~ 21 ' 21 

V is the relative velocity of the contact points defined in eq. ( |lq ) . 

We introduce relative coordinates r 12 = ri — r 2 and r — r x and the 
variables 

Z := ri2 • u± , 

Ml • M 2 J. 

a := r i2 • «i =w T 'i2 ' M i = ~ s i2 - 

y/\ - (ui ■ u 2 y 
1 x 

r i2 • «! = S 2 1 • 



^/l - («1 -M 2 ) 2 



The Jacobian of the transformation is given by \J\ — [u\ ■ U2) 2 . We remark 
that 4z Ir^l = V ■ M^sign(ri2 • u ± ) and we find again the relative velocity 
of the contact points V — ^12- — au\ — bii 2 given in the new coordinates. 
Integration over r gives V and integration over z gives the sum of two &— 
functions 0(±V ■ u±). This reflects the fact that if one particle touches the 
other from 'above' the sign of the relative velocity of the contact point has 
to be negative, if the particle touches from 'below' the velocity has to be 
positive. Next one introduces relative and center of mass momenta as well as 
dimcnsionless variables: 

1 1 r 

X:=-===(pi-j>2), !■= ^^^ =(pi+p 2 ), 

\l iml trans V^ mJ trans 

P6, - PA, 

Pd z ■= —?=== , VAi ■ = 



y/IT Iot ' \/IT mt sin Oi 

The integration over 7 can be done and the result is proportional to 



y / da db d(f>i sin^iii^i d(f>2 sin^2^2 d% dpg 1 dp t p 1 dpg 2 dp^ 

LI *-' 

v/l-(Mi-M 2 ) 2 cxp[--(x 2 + p* l + p\ 2 +pl +pl 2 )] 



p=±l 



V-ui 



0[p 



Vu, 



0(\a\-L/2)0(\b\-L/2)AE 12 , (80) 



all expressed in new variables and u by eq. (|74[), e.g. 



V = \/2T trans /mx - ay / T mt /I(p 9l e 9l + p Al e Al ) 



by / T mt /I{pe 2 eg 2 + p^e^) . (81) 
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We want to perform the remaining Gaussian integrals, but we have expressed 
different terms either in (u^, u ± ) defined according to eq. jl^ ) with i = 1,2 
or in (eg^e^J. It is useful to note that (u^,u ) and (e$ i ,e,j >i ) are two 
different orthonormal basis of the plane perpendicular to Ui, so that we can 
make a orthogonal coordinate transformation from one system to the other. 
The variables pg i and p ( j >1 are now standard normally distributed and after 
a orthogonal coordinate transformation the new coordinates will again be 
standard normally distributed. This means we can equivalently write (pe i ee i + 
P<t>i e <t>i) as ( v i u i~ +u>iU' L ) with standard normally distributed variables Vi and 
Wi. With this definition of vi and Wi we are able to evaluate for example terms 
of the form {pQ 1 ee 1 +p$ 1 e$ 1 ) -u± = (viu^+wiu±) -u± = »i, where we used 
uf- ■ u = u^ ■ u = 0. We can integrate freely over v\ and V2 and the two 
components of x perpendicular to u±. We denote with dQi :— d<pi sm(9i)d9i 
and the intermediate result reads 



^ 2V (Airy (2tt)(3/2) / ^ db dS dQxdn * exp (~ 2^ ) ^ 1 ~ (Ml ' U ^ 
G-s\0( P G-s) 



p=±i 



2T tI l + e„ 1 

1 s n 2 J- + g. + |i 



2/ T 2/ 




(82) 



2/ 



We introduced the vectors s := (si,S2,S3) := {x ' "1,101,11)2) and G 
^,-a\/%t,-&i 




We can now perform the integral over s. We sketch here only how this is 
done. We want to integrate 

[ ds exp(-is 2 )6)(±G • s)\G ■ s\(G ■ s)s 1 . (83) 

Let (ei, e2, 63) be the original coordinate system and we define a coordinate 
system (e x , e y , e z ) in which the z-axis is parallel to G and we decompose s 
in this coordinate system s = (s x , s y , s z ). Then eq. ( J83] ) reads 

/ ds x ds y ds z exp(--(sl + s 2 y + s 2 z ))0{±s z ) 

\G\\s z \\G\s z [(s x e x + s y e y + s z e z ) ■ e^ . (84) 

Only the term, which is proportional to s z e z , contributes and the Gaussian 
integral can easily be performed. Using that \G\e z = G we write \G\e z ■ e\ = 
G ■ e\= G\ and we end up with the result 4w\G\Gi. Only the integrals over 
Q\ and i?2 have to be done with standard techniques. All other integrals are 
performed similarly and the results are quoted in the main text. 
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